1. Stats de base, modèles linéaires, GLM

Distribution statistique

Décrit les probabilités d’apparition des valeurs d’une distribution donnée.

Distribution normale

On parle aussi souvent d’une distribution Gaussienne.


\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \]


\(\mu\) et \(\sigma\) sont des paramètres, respectivement la moyenne et l’écart-type.

Forme de la distribution normale

x <- seq(-5, 15, by = 0.01)
mu <- 5
sigma <- 2
y <- (1/(sigma * sqrt(2 * pi)) * exp((-1/2) * ((x - mu)/sigma)^2))
plot(x, y, type = "l")

Simuler une distribution normale

On peut utiliser la fonction rnorm pour simuler des valeurs provenant d’une distribution normale.

set.seed(123)
n <- 1000
x <- rnorm(n, mu, sigma)
hist(x, breaks = 50)

Interprétation

v <- seq(-5, 15, by = 0.01)
hist(x, breaks = 50, freq = FALSE, xlim = range(v))
curve(dnorm(x, mu, sigma), add = TRUE)
y <- dnorm(v, mu, sigma)
polygon(c(v[v >= 6 & v <= 8], 8, 6), c(y[v >= 6 & v <= 8], 0, 0), col = alpha("red", 0.5), border = NA)

L’aire sous la courbe est de 1 et la proportion de la surface en rouge représente la probabilité d’obtenir une valeur entre 6 et 8 avec une moyenne de \(\mu = 5\) et un écart-type de \(\sigma = 2\). On peut estimer cette probabilité avec la fonction pnorm intégrée à R. Ce type de fonction est intégré pour plusieurs distribution statistique (voir ?Distributions).

Probabilité d’obtenir une valeur?
pnorm(5, mean = 5, sd = 2)  # probabilité d'obtenir une valeur <= 5
## [1] 0.5
pnorm(-2, mean = 5, sd = 2)  # probabilité d'obtenir une valeur <= -2
## [1] 0.0002326291
pnorm(4, mean = 5, sd = 2)  # probabilité d'obtenir une valeur <= 4
## [1] 0.3085375

Probabilité d’obtenir une valeur entre 6 et 8

pnorm(8, 5, 2) - pnorm(6, 5, 2)
## [1] 0.2417303

Distribution d’une variable


\[ X \sim \mathcal{N}(\mu,\,\sigma^{2})\,\]


\(X\) suit une distribution normale avec une moyenne \(\mu\) et une variance \(\sigma^{2}\)

Erreur-type


\[ se = \frac{sd}{\sqrt{n}} \]


ou


\[ \sigma_{\bar{x}} = \frac{\sigma_{x}}{\sqrt{n}} \]


n <- 1e+05
pop <- rnorm(n, 100, 20)  # population fictive de n individus avec une moyenne de 100 et écart-type de 20
ech <- sample(pop, 30)  # échantillon aléatoire de 30 individus
mean(ech)  # moyenne
## [1] 100.5529
sd(ech)/sqrt(30)  # erreur-type
## [1] 3.27161

Simulation

Ici, on refait cet échantillonnage nreps fois et on calcule l’écart-type des différentes moyennes obtenues.

nreps <- 1000
ech1000 <- sapply(1:nreps, function(i) {
    s <- sample(pop, 30)
    mean(s)
})
hist(ech1000)

Écart-type des valeurs simulées

sd(ech1000)
## [1] 3.734432

Comparons ceci à l’erreur-type calculé à partir de notre seul échantillon:

sd(ech)/sqrt(30)  # erreur-type
## [1] 3.27161
Interprétation de l’erreur-type

L’écart-type de ces moyennes correspond à l’erreur-type estimée à partir d’un seul échantillon. En d’autres mots, l’erreur-type représente:

  • l’écart-type si des moyennes obtenues si on refaisait plusieurs fois le même échantillonnage

  • une mesure de précision dans l’estimation de la moyenne


Le concept d’erreur-type s’applique aussi à d’autres paramètres estimés et non seulement à la moyenne.


L’erreur-type d’une statistique ou d’un paramètre est une estimation de l’écart-type de sa distribution d’échantillonnage.

Intervalle de confiance

Reprenons notre population fictive pop qui a une moyenne de 100 et un écart-type de 20. On se rappelle qu’un intervalle de confiance pour une moyenne provenant d’un échantillon distribué normalement (ou presque) se calcule de la façon suivante:


\[\bar{x} ± t_{dl,\alpha/2}\sigma_{\bar{x}}\]


On peut calculer cet intervalle de confiance avec R:

mean(ech) + qt(0.025, df = length(ech) - 1, lower.tail = FALSE) * (sd(ech)/sqrt(length(ech))) * c(-1, 1)
## [1]  93.86168 107.24407

Simulation

Maintenant, reprenons note population fictive et calculons à chaque fois cet intervalle de confiance. Ensuite, calculons quelle est la proportion des intervalles qui contiennent la véritable moyenne de 100.

nreps <- 1000
ech1000 <- lapply(1:nreps, function(i) {
    s <- sample(pop, 30)
    c(mean(s), mean(s) + qt(0.025, df = length(s) - 1, lower.tail = FALSE) * (sd(s)/sqrt(length(s))) * c(-1, 1))
})
ci <- as.data.frame(do.call("rbind", ech1000))
names(ci) <- c("mean", "lowerCI", "upperCI")

sum(100 >= ci$lowerCI & 100 <= ci$upperCI)/nreps
## [1] 0.952

Interprétation de l’intervalle de confiance

Un intervalle de confiance à 95% veut dire que 95% des intervalles de confiance qu’on obtiendrait si on refaisait le même échantillonnage contiendront la véritable moyenne de la population.

Plus généralement, on devrait parler d’un paramètre (moyenne, variance, pente, etc.) plutôt que d’une moyenne.

On entend souvent à tort que cela veut dire qu’il y a 95% de chances que le paramètre soit à l’intérieur des bornes. Or le paramètre est considéré comme étant fixe et la probabilité porte sur le fait que l’intervalle contienne la véritable valeur du paramètre ou pas.

Tout comme l’erreur-type, l’intervalle de confiance est une mesure de précision dans l’estimation d’un paramètre.

Valeur de p

Un exemple avec le test de \(t\) et une population fictive.

set.seed(1)
n <- 10000
pop1 <- rnorm(n, 20, 3)  # population fictive de n individus avec une moyenne de 20
pop2 <- rnorm(n, 21, 3)  # population fictive de n individus avec une moyenne de 21

ech1 <- sample(pop1, 10)  # échantillon de 10 individus de chaque population
ech2 <- sample(pop2, 10)

test <- t.test(ech1, ech2)
test
## 
##  Welch Two Sample t-test
## 
## data:  ech1 and ech2
## t = -1.189, df = 14.686, p-value = 0.2533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.030060  1.147318
## sample estimates:
## mean of x mean of y 
##  19.22124  20.66261

Interprétation de la Valeur de p

plot(0, 0, xlim = c(-4, 4), ylim = c(0, 0.5), type = "n")
curve(dt(x, test$parameter), add = TRUE)
v <- seq(-5, 15, by = 0.01)
y <- dt(v, test$parameter)
polygon(c(v[v <= test$statistic], test$statistic, min(v)), c(y[v <= test$statistic], 0, 0), col = alpha("red", 0.5), border = NA)
polygon(c(v[v >= abs(test$statistic)], max(v), abs(test$statistic)), c(y[v >= abs(test$statistic)], 0, 0), col = alpha("red", 0.5), border = NA)
abline(v = rep(test$statistic, 2) * c(1, -1), lwd = 2, lty = 3, col = "red")

Simulations

Reprenons le test de \(t\) précédent, mais cette fois, prenons nos échantillons dans la même population (pop1) à chaque fois.

nreps <- 1000

pvalue <- sapply(1:nreps, function(i) {
    ech1 <- sample(pop1, 5)
    ech2 <- sample(pop1, 5)
    test <- t.test(ech1, ech2)
    test$p.value  #statistic
})

hist(pvalue, breaks = 50)
abline(v = 0.05, col = "red", lwd = 2, lty = 3)

sum(pvalue <= 0.05)/length(pvalue)
## [1] 0.051

Modèles linéaires simples


On a ici la formulation classique ou habituelle d’un modèle linéaire simple.


\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n} + \epsilon \]

\[ \epsilon \sim \mathcal{N}(0,\,\sigma^{2})\ \]

Autre formulation


\[ \mu = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n}\]

\[ y \sim \mathcal{N}(\mu,\,\sigma^{2})\ \]



\(y\) = observations

\(\mu\) = prédicteur linéaire

On a donc un modèle représentant le lien entre nos variables et on cherche à estimer les paramètres de ce modèle à l’aide de données.

Suppositions de bases

Quelles sont les suppositions de bases d’un modèle linéaire tel que celui présenté?

  • Normalité des résidus?

  • Homoscédasticité de la variance (la variance ne dépend pas de ou ne varie pas avec \(x_{n}\))

  • Linéarité des relations entre \(x_{n}\) et \(y\)?

  • Corrélation entre les \(x_{n}\)?

  • Indépendence entre les observations?

Simulations

  • Une des meilleures façons de s’assurer qu’on comprend comment fonctionne un modèle

  • Permet d’étudier le comportement d’un modèle avec des valeurs connues

Simulons un modèle

# nb d'observations
n <- 75

# valeurs fictives des variables explicatives selon une distribution uniforme
x1 <- runif(n, 0, 100)
x2 <- runif(n, 0, 10)
x3 <- runif(n, 0, 1)

# valeurs des paramètres beta
b0 <- 100
b1 <- 0.2
b2 <- 1
b3 <- -10

# prédicteur linéaire
mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3

# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)

# données
d <- data.frame(y, x1, x2, x3)

Fonction lm

La fonction lm permet de faire des modèles linéaires simples, ce qui inclue les régressions simples, les ANOVAs, ANCOVAs, régression multiples, etc.

Remarquez l’utilisation de l’argument data où la fonction ira chercher les différents éléments spécifiés dans la formule.

m <- lm(y ~ x1 + x2 + x3, data = d)

Coefficients

On peut appliquer plusieurs fonctions pour extraire les éléments d’un modèle, en particulier un modèle effectué avec la fonction lm.

coef(m)
## (Intercept)          x1          x2          x3 
##  96.9666068   0.1882304   1.0717263  -4.0315717

Sommaire

La fonction summary est la plus utile pour extraire l’information d’un modèle.

summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.518  -2.927   0.206   2.739  12.382 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 96.96661    1.80048  53.856  < 2e-16 ***
## x1           0.18823    0.01909   9.858 6.23e-15 ***
## x2           1.07173    0.18827   5.693 2.62e-07 ***
## x3          -4.03157    1.92616  -2.093   0.0399 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.02 on 71 degrees of freedom
## Multiple R-squared:  0.6444, Adjusted R-squared:  0.6294 
## F-statistic: 42.89 on 3 and 71 DF,  p-value: 6.301e-16

Vérification des suppositions de bases

On peut utiliser la fonction plot qui lorsqu’elle est appliquée sur un objet de classe lm produira des graphiques permettant de vérifier certaines suppositions de bases ou conditions d’application.

par(mfrow = c(2, 2))
plot(m)

Fonctions accessoires

On peut également créer facilement certains de ces graphiques avec d’autres fonctions.

par(mfrow = c(1, 3))
plot(resid(m), fitted(m))
hist(resid(m))
qqnorm(resid(m))
qqline(resid(m))

Facteurs

# simulation d'un facteur
d$fac <- factor(sample(c("A", "B", "C"), nrow(d), replace = TRUE))
d$y <- d$y + 2 * as.integer(d$fac) - 1

# modèle modifié avec le facteur
m <- lm(y ~ x1 + x2 + x3 + fac, data = d)
summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + fac, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.3549  -3.0036   0.3166   2.6305  13.2288 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  97.2099     1.9870  48.924  < 2e-16 ***
## x1            0.1861     0.0190   9.794 1.09e-14 ***
## x2            1.0464     0.1854   5.645 3.39e-07 ***
## x3           -3.8057     1.8956  -2.008   0.0486 *  
## facB          1.9366     1.4416   1.343   0.1836    
## facC          6.5096     1.3744   4.736 1.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.932 on 69 degrees of freedom
## Multiple R-squared:  0.696,  Adjusted R-squared:  0.674 
## F-statistic:  31.6 on 5 and 69 DF,  p-value: < 2.2e-16

Comprendre son modèle avec des graphiques

La fonction predict sert à calculer la valeur prédite par le modèle pour chaque observation. Par exemple, ceci permet de mettre en relation les valeurs observées et les valeurs prédites. La fonction predict dans ce cas calcule la valeur prédite par le modèle pour chaque observation dans la base de donnée.

plot(d$y, predict(m))

Soumettre de nouvelles valeurs à predict

Plus souvent, on veut illustrer l’effet de nos variables et pour cela, il faut soumettre des valeurs à la fonction predict. Cette fonction est très utile lorsque l’on veut un maximum de contrôle sur les valeurs à prédire.

x <- seq(min(d$x1), max(d$x1), length.out = 50)
nd <- data.frame(x1 = x, x2 = mean(d$x2), x3 = mean(d$x3), fac = "B")  # on fixe les autres variables à leur valeur moyenne
p <- predict(m, newdata = nd)
plot(d$x1, d$y)  # observations
lines(x, p)  # valeurs prédites

Package visreg

Le package visreg permet d’illustrer facilement les prédictions d’un modèle (ou les effets marginaux des différentes variables explicatives). Il utilise en arrière-plan la fonction predict.

library(visreg)
par(mfrow = c(2, 2))
visreg(m)

Package ggeffects

Le package ggeffects permet également d’illustrer les prédictions d’un modèle, mais ce package se base sur l’utilisation du package ggplot2 pour construire les graphiques. Le principe demeure le même: en arrière-plan la fonction predict est utilisée. On peut également combiner les différents graphiques générés en utilisant le package patchwork et sa fonction wrap_plots.

library(ggeffects)
library(patchwork)
g <- ggpredict(m)
p <- plot(g, add.data = TRUE)
wrap_plots(p)

Interactions

Dans vos mots, qu’est-ce qu’une interaction?


Un exemple

# modèle modifié avec le facteur
m <- lm(y ~ x1 + x2 + x3 + x1 * fac, data = d)
summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x1 * fac, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0406  -1.9485   0.1316   2.1335  14.2374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 98.50371    2.15978  45.608  < 2e-16 ***
## x1           0.14204    0.03151   4.508  2.7e-05 ***
## x2           1.15443    0.18388   6.278  2.9e-08 ***
## x3          -2.44076    1.86546  -1.308  0.19521    
## facB         0.66407    2.63788   0.252  0.80201    
## facC        -1.01455    3.01397  -0.337  0.73746    
## x1:facB      0.01671    0.04332   0.386  0.70090    
## x1:facC      0.13612    0.04871   2.794  0.00678 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.692 on 67 degrees of freedom
## Multiple R-squared:  0.7329, Adjusted R-squared:  0.705 
## F-statistic: 26.26 on 7 and 67 DF,  p-value: < 2.2e-16

Illustration graphique

visreg(m, "x1", by = "fac", overlay = TRUE)


On a une interaction significative, mais est-elle importante d’un point de vue biologique? Peut-être pas… Toujours distinguer entre la signification statistique et la signification biologique.

Une valeur de p significative c’est bien beau, mais ce qui importe vraiment c’est la taille de l’effet.

Interactions avec ggeffects
g <- ggpredict(m, terms = c("x1", "fac"), message = FALSE)
plot(g, add.data = TRUE, jitter = FALSE)

Interprétion

  • Beaucoup plus facile avec un graphique qu’en regardant une table de coefficients.

  • En présence d’une interaction (significative), les effets simples sont plus difficilement interprétables et ne peuvent être interprétés sans faire référence à l’interaction (à moins d’ajustements particuliers voir Schielzeth 2010).

  • Si pour deux variables impliquées dans une interaction les effets simples ne sont pas significatifs, mais que leur interaction l’est, il faut considérer que ces deux variables ont un effet même si les coefficients associés aux effets simples ne sont pas significatifs.

Types d’erreurs


  • Erreur de type I: rejeter l’hypothèse nulle, alors qu’elle est vraie (faux positif).


  • Erreur de type II: ne pas rejeter l’hypothèse nulle, alors qu’elle est fausse (faux négatif).

Simulations

Simulons un modèle dans lequel les variables n’ont aucun effet.

Que se passera-t-il si on répète plusieurs fois ce scénario?

À quoi ressembleront nos valeurs de p?

Créons d’abord une fonction sim_model pour générer un modèle linéaire fictif.

sim_model <- function(b0, b1, b2, b3) {
    n <- 200
    x1 <- runif(n, 0, 100)
    x2 <- runif(n, 0, 10)
    x3 <- runif(n, 0, 1)
    mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3
    y <- rnorm(n, mu, 5)
    d <- data.frame(y, x1, x2, x3)
    m <- lm(y ~ x1 + x2 + x3, data = d)
    m
}

Simulons

Ensuite, simulons nsims modèles fictifs et emmagasinons ces modèles dans une liste nommée list_models.

nsims <- 500

list_models <- lapply(1:nsims, function(i) {
    sim_model(b0 = 100, b1 = 0, b2 = 0, b3 = 0)
})


Si on illustre les valeurs de p obtenues pour chacun des coefficients associés aux variables \(x_{n}\) à l’aide d’un histogramme, à quoi ressemblera cet histogramme?

Valeurs de p obtenus

p <- lapply(list_models, function(i) {
    summary(i)$coef[2:4, 4]
})

hist(unlist(p), breaks = seq(0, 1, by = 0.05), xlab = "Valeurs de p")

Conclusions
  • Dans une certaine proportion des cas ( ~ 5% ), on conclut à un effet, alors qu’il n’y en a pas pas! C’est le fameux seuil \(\alpha\) de 0.05.

  • Sous l’hypothèse nulle, les valeurs de p sont distribuées uniformément entre 0 et 1.

Modèles linéaires standards

\[y \sim \mathcal{N}(\mu, \sigma^2)\]

\[\mu=\beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_k X_k\]

Formulation matricielle

\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{I}\sigma^2)\]


\(\beta\) est un vecteur de coefficients.

\(\mathbf{X}\) est une matrice avec les valeurs des variables explicatives.

\(\mathbf{I}\) est une matrice identité qui a pour effet d’assigner une variance égale \(\sigma^2\) à toutes les observations.

Modèles linéaires généraux

\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{\Sigma})\]

\(\mathbf{\Sigma}\) est une matrice de variance-covariance qui permet beaucoup plus de flexibilité en terme de variance et de dépendance entre les observations. En général, des paramètres seront associés à la matrice \(\mathbf{\Sigma}\) pour permettre différentes situations. Dans le cas d’un modèle linéaire standard, on a \(\mathbf{\Sigma} = \mathbf{I}\sigma^2\).


Ex.: variance différente pour les niveaux d’un facteur, variance augmentant avec certaines valeurs de \(x\), autocorrélation temporelle ou spatiale, etc.


Ces différentes situations peuvent être obtenues entre autres avec les packages nlme (fonctions gls ou lme) ou glmmTMB.

Modèles linéaires généralisés (GLM)

\[y \sim \mathcal{ExpFam}(\theta, ...)\]

\[\mathcal{g}(\mu)=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]

Structure d’un GLM

\[y \sim \mathcal{ExpFam}(\theta, ...)\]


\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]


\[ g(E(y)) = g(\mu) = \eta \]

Composante systématique

\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]

Le prédicteur linéaire \(\eta\) permet de relier les variables explicatives à la variable réponse.

Composante aléatoire

\[ y \sim \mathcal{ExpFam}(\theta) \]


Les observations \(y\) suivent une distribution particulière appartenant à la famille des distributions exponentielles avec des paramètres \(\theta\). La distributions normale, de Poisson, Binomiale et Gamma font partie de la famille des distributions exponentielles.

Fonction de lien

\[ g(E(y)) = g(\mu) = \eta \]


La fonction de lien (link function) est une transformation qui permet de relier la réponse attendue \(E(y)\) (i.e. la moyenne) au prédicteur linéaire.


C’est ce qui permet de relier les variables explicatives à la variable réponse. Pour chaque distribution qui pourrait être utilisée dans le cadre d’un GLM, une fonction de lien dite canonique est utilisée par défaut, mais il existe plusieurs fonctions de liens possibles.

GLM Poisson

\[ y \sim \mathcal{Pois}(\lambda) \]


\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]


On a donc des observations \(y\) qui suivent une distribution de Poisson avec un paramètre \(\lambda\). Une fonction de lien \(log\) permet de relier linéairement les variables explicatives à la valeur attendue \(E(y) = \lambda\) de la variable réponse.

Distribution de Poisson

Cette distribution est souvent utilisée pour décrire le nombre de fois qu’un événement survient dans un intervalle de temps ou d’espace.

\[f(k,\lambda) = Pr(y = k) = \frac{\lambda^{k}e^{-\lambda}}{k!}\]

Cette fonction permet de calculer la probabilité d’observer un événement \(k\) fois avec une distribution de Poisson ayant un paramètre \(\lambda\). La moyenne de cette distribution est égale à \(\lambda\).

Cette distribution a la propriété que \(\lambda = E(y) = Var(y)\). Autrement dit, la moyenne et la variance sont égales et correspondent au paramètre \(\lambda\) de la distribution.

D’où provient-elle?

Simulons une distribution de Poisson en distribuant aléatoirement 200 observations dans un carré 10 x 10. Ensuite, divisons ce carré en 100 cellules.

library(sf)

n <- 200
x <- runif(n, 0, 10)
y <- runif(n, 0, 10)

sfc <- st_sfc(st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 0)))))
grid <- st_make_grid(sfc, cellsize = 1)
plot(grid, axes = TRUE)
points(x, y)

Nb d’observations par cellule

Ensuite, calculons le nombre d’observations par cellule et illustrons la distribution du nombre d’observations pour chaque cellule. En moyenne, on devrait retrouver 2 observations par cellule (200 obs. / 100 cell. = 2).

o <- st_intersects(grid, st_as_sf(data.frame(x, y), coords = c("x", "y")))
co <- sapply(o, length)
hist(co, breaks = seq(-0.5, max(co) + 0.5, by = 1), freq = FALSE, main = "", xlab = "Nb d'observations")
lines(0:max(co), dpois(0:max(co), lambda = 2), type = "b", cex = 2)
legend("topright", lty = c(NA, 1), legend = c("Observée", "Théorique"), bty = "n", cex = 1, pch = c(22, 21), pt.bg = c("gray", "white"))

Ce serait la même chose si on simulait des observations aléatoires dans le temps et si on comptait le nb d’observations par intervalles de temps.

Simulons un modèle

Simulons un modèle à partir d’une distribution de Poisson.

n <- 100
x <- runif(n, 0, 10)
lambda <- exp(-2 + 0.4 * x)
y <- rpois(n, lambda)
d <- data.frame(y, x)

m <- glm(y ~ x, family = "poisson", data = d)
summary(m)
## 
## Call:
## glm(formula = y ~ x, family = "poisson", data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.46287    0.30996  -7.946 1.93e-15 ***
## x            0.45294    0.03784  11.971  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 289.290  on 99  degrees of freedom
## Residual deviance:  85.894  on 98  degrees of freedom
## AIC: 245.83
## 
## Number of Fisher Scoring iterations: 5

Coefficients

Les coefficients estimés sont très similaires à ce qui a été utilisé pour générer le modèle (-2 et 0.4).

coef(m)
## (Intercept)           x 
##  -2.4628709   0.4529379

Prédictions

Illustrons l’effet de la variable explicative.

visreg(m)

Sous quelle échelle?

visreg(m, scale = "response")

Échelle du prédicteur linéaire ou de la réponse

par(mfrow = c(1, 2))
visreg(m)
visreg(m, scale = "response")

Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction log permet de transformer les comptes allant de 0 à l’\(+\infty\) en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle log, les valeurs entre 0 et 1 sont négatives et les valeurs supérieures à 1 sont positives.

\[ log(1) = 0 \] \[ log(0.1) = -2.302586 \] \[ log(5) = 1.609438 \]

ggeffects

g <- ggpredict(m, terms = c("x [n=50]"))
plot(g, add.data = TRUE, jitter = FALSE)

predict

v <- seq(min(x), max(x), length.out = 50)
plot(y ~ x, data = d)
p <- predict(m, data.frame(x = v))
lines(v, p)

Argument type

Par défaut, l’argument type = "link" et il faut spécifier type = "response" pour obtenir les prédictions sous l’échelle de la réponse.

v <- seq(min(x), max(x), length.out = 50)
plot(y ~ x, data = d)
p <- predict(m, data.frame(x = v), type = "response")
lines(v, p)

Suppositions de base

On pourrait être tenté d’utiliser la fonction plot pour vérifier les conditions d’applications.

par(mfrow = c(2, 2), mar = c(4, 4, 2, 2))
plot(m)

Quelles sont les suppositions de bases?

\[ y \sim \mathcal{Pois}(\lambda) \]

\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]

Certaines suppositions sont les mêmes que pour un modèle linéaire standard avec une distribution, mais d’autres diffèrent.

Vérification avec DHARMa

Une façon de vérifier les supposition de bases des modèles est d’utliser le modèle pour simuler des observations et ensuite comparer ces observations simulées aux observations réelles. Si les observations réelles ne ressemblent pas aux observations simulées à partir du modèle, c’est que le modèle ne représente pas bien ce qui se passe avec les données réelles et donc que qqch devrait être modifié dans le modèle.

Le package DHARMa se base sur ces simulations pour comparer les observations attendues par le modèle aux observations réelles. Plus spécifiquement, la méthode utilise les résidus quantiles pour comparer chaque observation aux observations attendues pour cette observations. Les résidus quantiles sont rapportés entre 0 et 1 et si le modèle est approprié pour les données, on devrait s’attendre à ce que les résidus soient distributés uniformément entre 0 et 1 pour toutes les observations.

library(DHARMa)

s <- simulateResiduals(m)
plot(s)

Si le modèle est appropré pour les données, on devrait voir ici une distribution uniforme des points et il ne devrait pas y avoir de patrons particuliers dans la distribution des points dans le graphique de droite. En particulier, les lignes de prédictions illustrées devraient être alignées sur les trois lignes représentant les quantiles 25, 50 et 75%. Dans le graphique quantile-quantile de gauche, les points devraient être alignés sur la droite.

Autres vérifications

Les problèmes de sous- ou sur- dispersion et de surabondance de 0 (zéro-inflation) sont très fréquents avec les GLM. Ces deux fonctions de DHARMa permettent d’évaluer si les observations réelles diffèrent des observations simulées en terme de présence de 0 ou de dispersion. En rouge, ce sont les observations réelles et en gris les observations simulées. Si les observations réelles sont bien en dehors des observations simulées, c’est qu’il y a un problème avec le modèle et il faut un ajustment.

par(mfrow = c(1, 2))
testDispersion(s)
testZeroInflation(s)

GLM binomial

\[ y \sim \mathcal{Binom}(n,p) \]

\[ logit(p) = log\left(\frac{p}{1-p}\right) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]

Régression logistique

\[ y \sim \mathcal{Binom}(1,p) = \mathcal{Bernoulli}(p)\]

La distribution de Bernoulli est un cas particulier de la distribution binomiale avec un tirage (e.g. un individu est vivant ou mort, un individu est mâle ou femelle, etc.)

Distribution binomiale

\[f(k,n,p) = Pr(y = k) = \binom{n}{k}p^k(1-p)^{n-k}\]

\[\binom{n}{k}=\frac{n!}{k!(n-k)!}\]


Ce qui représente le nombre de façons de prendre \(k\) éléments parmi \(n\) éléments.


Cette fonction permet de calculer la probabilité d’observer \(k\) succès si \(n\) tirages sont faits et que la probabilité de succès est de \(p\)


Cette distribution a la propriété que \(E(y) = np\) et \(Var(y) = np(1-p)\) , autrement dit, la moyenne et la variance sont entièrement déterminées par la probabilité de succès et par le nombre de tirages. On a donc une autre distribution assez restrictive.

Exemple

Supposons qu’on tire pile (0) ou face (1) 10 fois et qu’on calcule le nombre de face et qu’on répète cette exercice 500 fois. Quelle sera la distribution du nombre de faces obtenus?

n <- 500
co <- sapply(1:n, function(i) {
    s <- sample(0:1, 10, prob = c(0.5, 0.5), replace = TRUE)
    sum(s)  # on additione les faces (1)
})

hist(co, breaks = seq(-0.5, max(co) + 0.5, by = 1), freq = FALSE, main = "", xlab = "Nb de faces")
lines(0:max(co), dbinom(0:max(co), prob = 0.5, size = 10), type = "b", cex = 2)
legend("topright", lty = c(NA, 1), legend = c("Observée", "Théorique"), bty = "n", cex = 1, pch = c(22, 21), pt.bg = c("gray", "white"))

Simulons un modèle

logit <- function(p) {
    log(p/(1 - p))
}

inv.logit <- function(x) {
    exp(x)/(1 + exp(x))
}

n <- 200
x <- runif(n, 0, 10)
z <- -2 + 0.5 * x
y <- rbinom(n, size = 1, prob = inv.logit(z))

d <- data.frame(y, x)

m <- glm(y ~ x, family = "binomial", data = d)
summary(m)
## 
## Call:
## glm(formula = y ~ x, family = "binomial", data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.68335    0.43063  -6.231 4.63e-10 ***
## x            0.63864    0.08853   7.214 5.45e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 272.74  on 199  degrees of freedom
## Residual deviance: 185.42  on 198  degrees of freedom
## AIC: 189.42
## 
## Number of Fisher Scoring iterations: 5

Visualisation

Avec visreg

par(mfrow = c(1, 2))
visreg(m)
visreg(m, scale = "response")

Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction logit permet de transformer les probabilités allant de 0 à 1 en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle logit, une valeur de 0 corrrespond à une probabilité de 0.5.

\[ log\left(\frac{p}{1-p}\right) = log\left(\frac{0.5}{1-0.5}\right) = log(1) = 0 \]

Visualisation

Avec ggeffects

g <- ggpredict(m, terms = c("x [n=50]"))
plot(g, add.data = TRUE)

Suppositions de base

par(mfrow = c(2, 2))
s <- simulateResiduals(m)
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Problèmes

Les problèmes les plus fréquents sont la surdispersion et la surabondance de zéros (zero-inflation)

Surdispersion

C’est lorsque les données sont plus variables que ce qui est attendu par le modèle. Les distribution de Poisson et binomiale sont très restrictives par rapport à leur variance et les problèmes de surdispersion sont très fréquents.


En général, une surdispersion ignorée fera en sorte que l’incertitude sera sous-estimée (trop de confiance dans nos résultats).


Plusieurs causes possibles (variables explicatives manquantes, phénomènes d’aggrégation, hétérogénéité, etc.)


Solutions:

Poisson -> Négative binomiale (glm.nb, glmmTMB, brms, etc.)

Binomiale -> Béta binomiale (glmmTMB, brms, etc.)

Obervation-Level Random Effects (OLRE) (voir Harrison 2015 et Harrison 2014).

Exemple

Voici un exemple de GLM de Poisson avec de la surdispersion. Pour créer cette surdispersion, j’utilise une distribution négative binomiale pour générer les observations au lieu d’une distribution de Poisson.

n <- 100
x <- runif(n, 0, 10)
lambda <- exp(-2 + 0.4 * x)
y <- rnbinom(n, mu = lambda, size = 0.5)  # 
d <- data.frame(y, x)

m <- glm(y ~ x, family = "poisson", data = d)
g <- ggpredict(m, terms = "x")
plot(g, add.data = TRUE)

Vérification avec DHARMa

Les résidus ne sont clairement pas distributés uniformément sur le second graphique et la surdispersion est évidente sur le 3e graphique.

s <- simulateResiduals(m)
par(mfrow = c(1, 4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Correction avec glmmTMB

library(glmmTMB)

m <- glmmTMB(y ~ x, family = "nbinom1", data = d)

s <- simulateResiduals(m)

par(mfrow = c(1, 4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Sousdisperion

C’est lorsque les données sont moins variables que ce qui est attendu par le modèle. Ceci est beaucoup plus rare et probablement moins dommageable (e.g. taille de couvée/portée chez les oiseaux/mammifères voir Brooks et al. 2019).


En général, une sous-dispersion ignorée fera en sorte que l’incertitude sera surestimée (trop conservateurs dans nos résultats) ce qui est un moindre mal.


Solutions:

Poisson -> glmmTMB, family = "compois"

Zéro-inflation

Surabondance de 0 dans les observations par rapport à ce qui est attendu par le modèle. Il est possible aussi d’avoir une sous-abondance de 0 et/ou encore une surabondance de 1, mais ce dernier cas est un peu plus rare.


Solutions:

modèle zero-inflated (glmmTMB,pscl,brms, etc.)


À noter qu’une surdipersion peut être due à une surabondance de 0 et vice-versa. Ce sont deux problèmes qui peuvent interagir.

Extensions


Gamma    \(]0,\infty\)

Variable réponse continue, strictement positive et avec une asymétrie à droite.

Souvent une alternative à une transformation \(log\).

Ex.: concentration d’un composé, biomasse, etc.


Tweedie    \([0,\infty\)

Variable réponse continue positive et incluant des 0.

Ex.: biomasse capturée, précipitations, etc.

Voir Lecomte et al. 2013


Beta    \(]0,1[\)   ou   \([0,1]\)

Variable continue entre 0 et 1 exclusivement.

Pour proportion ou % ne provenant pas d’un ratio entre entiers.

Si 0 et 1 possibles, certaines transformations peuvent être utilisées.

Ex.: proportion d’une surface affectée, proportion de temps passé en alimentation, etc.

Voir Douma & Weedon 2019

2. Modèles mixtes/hiérarchiques



\[ y = \alpha + \beta x + \epsilon \]

\[ \epsilon \sim \text{Normale}(0,\sigma^{2}) \]

ou


\[ \mu = \alpha + \beta x \]

\[ y \sim \text{Normale}(\mu,\sigma^{2}) \]

Modèle linéaire hiérarchique simple



\[ y_{ij} = \alpha + \beta x_{ij} + \tau_{i} + \epsilon_{ij}\]

\[ \tau_{i} \sim \text{Normale}(0,\,\sigma_{\tau}^{2}) \]

\[ \epsilon_{ij} \sim \text{Normale}(0,\sigma_{\epsilon}^{2})\ \]



où il y a \(i\) groupes et \(j\) observations par groupe

Autre formulation



\[ \mu_{ij} = \alpha_{i} + \beta x_{ij}\]

\[ \alpha_{i} \sim \text{Normale}(\alpha,\,\sigma_{\alpha}^{2})\ \]

\[ y_{ij} \sim \text{Normale}(\mu_{ij},\,\sigma_{\epsilon}^{2})\ \]

Un modèle est dit hiérarchique lorsque certains de ses paramètres sont définis par une distribution de probabilité.

Simulons un exemple

set.seed(12345)
ng <- 20
nobs <- 10
n <- ng * nobs
id <- gl(ng, nobs)

# valeurs fictives des variables explicatives selon une distribution uniforme
x <- runif(n, 0, 100)

# valeurs des paramètres
a <- 100
b <- 1

# on génère l'effet aléatoire pour chaque groupe
rea <- rep(rnorm(ng, 0, 20), each = nobs)

# prédicteur linéaire
mu <- (a + rea) + b * x

# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)

# données
d <- data.frame(y, x)

Illustrons le modèle

library(lme4)
library(ggeffects)

m <- lmer(y ~ x + (1 | id), data = d)

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)

Illustrons les pentes par individus

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = sample(colors(), 20))

Retrouvons les paramètres

La variance de l’effet aléatoire était de 20² et la variance résiduelle était de 5².

summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | id)
##    Data: d
## 
## REML criterion at convergence: 1314.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50905 -0.67270 -0.04039  0.61012  2.58655 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 386.9    19.67   
##  Residual              25.3     5.03   
## Number of obs: 200, groups:  id, 20
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 101.35098    4.46512   22.70
## x             0.99676    0.01262   78.97
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.153

Récupérons ces effets

fixef(m)
## (Intercept)           x 
## 101.3509752   0.9967622
ranef(m)
## $id
##    (Intercept)
## 1    2.8823702
## 2  -22.7874952
## 3    5.4199340
## 4  -28.1373544
## 5    1.0787346
## 6  -10.6043879
## 7   -6.4207965
## 8   31.3275567
## 9  -10.5628148
## 10   4.4707573
## 11 -24.5843971
## 12 -26.8974001
## 13  24.9063468
## 14  23.4502635
## 15  -2.0989706
## 16 -12.2974161
## 17  -1.4157682
## 18   8.9803887
## 19  44.1496752
## 20  -0.8592262
## 
## with conditional variances for "id"
apply(ranef(m)$id, 2, sd)
## (Intercept) 
##    19.60505

Coefficients

coef(m)
## $id
##    (Intercept)         x
## 1    104.23335 0.9967622
## 2     78.56348 0.9967622
## 3    106.77091 0.9967622
## 4     73.21362 0.9967622
## 5    102.42971 0.9967622
## 6     90.74659 0.9967622
## 7     94.93018 0.9967622
## 8    132.67853 0.9967622
## 9     90.78816 0.9967622
## 10   105.82173 0.9967622
## 11    76.76658 0.9967622
## 12    74.45358 0.9967622
## 13   126.25732 0.9967622
## 14   124.80124 0.9967622
## 15    99.25200 0.9967622
## 16    89.05356 0.9967622
## 17    99.93521 0.9967622
## 18   110.33136 0.9967622
## 19   145.50065 0.9967622
## 20   100.49175 0.9967622
## 
## attr(,"class")
## [1] "coef.mer"

Complexifions le modèle

b
## [1] 1
reb <- rep(rnorm(ng, 0, 0.5), each = nobs)

mu <- (a + rea) + (b + reb) * x

y <- rnorm(n, mu, 5)

d <- data.frame(y, x)

Illustrons les pentes

m <- lmer(y ~ x + (x | id), data = d)

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)

Retrouvons les paramètres

La variance de l’effet aléatoire sur la pente était de (0.5)².

summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (x | id)
##    Data: d
## 
## REML criterion at convergence: 1370.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.19310 -0.58325  0.02874  0.58522  2.42595 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 388.2872 19.7050       
##           x             0.2105  0.4588  -0.27
##  Residual              22.0018  4.6906       
## Number of obs: 200, groups:  id, 20
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 100.4564     4.4699  22.474
## x             0.8959     0.1033   8.672
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.283
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00784532 (tol = 0.002, component 1)

Récupérons ces effets

fixef(m)
## (Intercept)           x 
## 100.4563555   0.8958772
ranef(m)
## $id
##    (Intercept)           x
## 1    -4.259743 -0.11207853
## 2   -24.555894  0.09388618
## 3     9.760813  0.08930651
## 4   -23.972706 -0.43249104
## 5     6.379230  0.06853379
## 6   -10.672292  0.92265805
## 7    -4.656865 -0.14388477
## 8    33.097624 -0.22946144
## 9    -7.593899  0.23044679
## 10   -2.310431  0.38486838
## 11  -32.815439 -0.47826203
## 12  -26.583571  0.53612804
## 13   22.768061  0.34741573
## 14   22.280452 -0.58389607
## 15    1.432798  0.75717249
## 16   -8.201353 -0.31288440
## 17    1.063482  0.10990165
## 18   10.101182  0.03586551
## 19   39.745709 -0.88012487
## 20   -1.007156 -0.40309995
## 
## with conditional variances for "id"
apply(ranef(m)$id, 2, sd)
## (Intercept)           x 
##  19.4359900   0.4557504
coef(m)
## $id
##    (Intercept)          x
## 1     96.19661 0.78379868
## 2     75.90046 0.98976339
## 3    110.21717 0.98518372
## 4     76.48365 0.46338617
## 5    106.83559 0.96441100
## 6     89.78406 1.81853525
## 7     95.79949 0.75199243
## 8    133.55398 0.66641577
## 9     92.86246 1.12632400
## 10    98.14592 1.28074559
## 11    67.64092 0.41761518
## 12    73.87278 1.43200525
## 13   123.22442 1.24329294
## 14   122.73681 0.31198114
## 15   101.88915 1.65304970
## 16    92.25500 0.58299281
## 17   101.51984 1.00577886
## 18   110.55754 0.93174272
## 19   140.20206 0.01575234
## 20    99.44920 0.49277726
## 
## attr(,"class")
## [1] "coef.mer"

Formulations

m <- lmer(y ~ x + (1 | id), data = d)

m <- lmer(y ~ x + (x | id), data = d)  # équivalent à (1 + x|id)

m <- lmer(y ~ x + (0 + x | id), data = d)

m <- lmer(y ~ x + (1 | id) + (0 + x | id), data = d)

Corrélations intercepts/pentes

m <- lmer(y ~ x + (1 | id) + (0 + x | id), data = d)


Le dernier est équivalent à

m <- lmer(y ~ x + (x || id), data = d)


et signifie qu’on assume que l’intercept et la pente ne sont pas corrélées ou qu’on ne veut pas estimer cette corrélation. En pratique, il n’y a pas énormément de raison de faire cela et par défaut le modèle va estimer cette corrélation.

Effets aléatoires croisés

Supposons que nous avons des individus observés plusieurs fois par année et ce sur plusieurs années.

m <- lmer(y ~ x + (1 | id) + (1 | year), data = d)

Effets aléatoires nichés

Supposons maintenant que les individus sont annuels et qu’ils ne sont observés que pour des années particulières.

m <- lmer(y ~ x + (1 | year/id), data = d)


Équivalent à:

m <- lmer(y ~ x + (1 | year) + (1 | year:id), data = d)

Quand utiliser les effets aléatoires?


  • Mesures répétées sur certaines unités ou individus

  • Les observations sont regroupées sous des structures ou ne sont pas indépendantes

  • Il y a de la réplication à plusieurs niveaux de la hiérarchie

  • On veut estimer la variance ou généraliser à une population d’unités d’échantillonnage

Effet spatial, temporel, phylogénétique


Autocorrélation spatiale, temporelle, phylogénétique, etc…

Régularisation (shrinkage)


Complete pooling

No pooling

Partial pooling

Un exemple

Effet du n

Problèmes, questionnements ?



La GLMM FAQ!


https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html

3. Vraisemblance et statistiques bayésiennes

Rappel des probabilités

Quelle est la probabilité d’obtenir 2 fois pile en tirant une pièce de monnaie 2 fois?

Pile ou face

On doit multiplier la probabilité de chaque événement pour obtenir la probabilité combinée de ces deux événements. Si \(P(pile) = 0.5\) (i.e. la probabilité d’obtenir pile), on a:


\[P(pile,pile) = P(pile) \times P(pile) = 0.5 \times 0.5 = 0.25\]

Avec R

On peut calculer ceci avec R et la fonction dbinom.

dbinom(0, size = 1, p = 0.5) * dbinom(0, size = 1, p = 0.5)
## [1] 0.25


On peut également appliquer la fonction dbinom à un vecteur de valeur. Ceci nous donne un vecteur de (densité de) probabilité et on peut multiplier ces valeurs entre elles avec la fonction prod.

prod(dbinom(c(0, 0), size = 1, p = 0.5))
## [1] 0.25

Vraisemblance

La vraisemblance (likelihood) représente la plausibilité de valeurs de paramètres en fonction de données observées. L’utilisation de la vraisemblance implique l’utilisation de distribution de probabilités pour représenter le lien entre les observations et les modèles.

Pour quantifier la vraisemblance, on peut utiliser une fonction de vraisemblance (likelihood function).
\[L(\theta|data) = P(data|\theta)\]
\(data\): données
\(\theta\): paramètres du modèle

\(L(\theta|data)\): fonction de vraisemblance
\(P(data|\theta)\): vraisemblance (likelihood)


Ici, \(L(\theta|data)\) est la fonction de vraisemblance et l’équation précédente nous dit que la vraisemblance des paramètres conditonnelle aux paramètres \(\theta\) est égale à la (densité de) probabilité des données conditionnelle aux paramètres \(\theta\).

Maximum de vraisemblance

L’idée derrière le maximum de vraisemblance (maximum likelihood) est de déterminer la valeur des paramètres qui rend l’observation des données la plus probable.



Une grande partie des paramètres que nous estimons sont estimés par la méthode du maximum de vraisemblance.
Certaines méthodes plus classiques comme la minimisation de la somme des carrés dans la régression standard sont équivalentes à la maximisation de la vraisemblance. La sélection de modèles par AIC repose sur l’utilisation du maximum de vraisemblance.

Exemple

Simulons n = 20 présences/absences avec une probabilité de p = 0.5. Ensuite, calculons la vraisemblance d’observer les données générées en fonction de plusieurs valeurs du paramètre p, qui ici est une probabilité.


set.seed(123)
p <- 0.5
y <- rbinom(20, 1, p)  # tirage aléatoire de 20 valeurs de 0/1 avec une probabilité de 0.5

probs <- seq(0, 1, by = 0.01)  # valeurs du paramètre qui seront évaluées
l <- sapply(probs, function(i) {
    prod(dbinom(y, size = 1, prob = i))  # likelihood
})


Pour calculer la vraisemblance, on utilise la fonction dbinom qui permet de calculer pour chaque valeur observée une (densité de) probabilité. En calculant le produit prod de toutes ces valeurs, on obtient la valeur de la fonction de vraisemblance pour l’observation de l’ensemble de ces données pour une valeur donnée du paramètre p. On calcule ensuite ces valeurs de vraisemblance pour un ensemble de valeur de p.

Profil de vraisemblance

Par la suite, on peut illustrer la vraisemblance d’obtenir les valeurs observées pour les différentes valeurs de p. Ceci permet d’illustrer le profile de vraisemblance (likelihood profile).

plot(probs, l, type = "l", ylab = "Likelihood", xlab = "Valeur du paramètre p")
abline(v = sum(y)/length(y), lty = 3)
abline(v = 0.5)
legend("topleft", lty = c(1, 3), legend = c("Probabilité réelle", "Proportion de 1 dans l'échantillon"), bty = "n")

probs[which.max(l)]
## [1] 0.55


Le maximum de cette courbe est le maximum de vraisemblance (maximum likelihood), i.e. la valeur de p pour laquelle la probabilité d’observer les données est la plus grande. Dans ce cas précis, cette valeur est la proportion de 1 dans les données. La plupart des paramètres que nous estimons sont obtenus avec cette méthode.

Avec un glm

Voici ce qu’on obtient avec un glm si on tente d’estimer cette probabilité. En principe, l’intercept représente cette probabilité, mais elle est sur l’échelle logit \(log(p/(1-p))\). Il faut donc retransformer le tout sous l’échelle d’une probabilité avec une fonction logit inverse.

m <- glm(y ~ 1, family = binomial)
summary(m)
## 
## Call:
## glm(formula = y ~ 1, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.2007     0.4495   0.446    0.655
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.526  on 19  degrees of freedom
## Residual deviance: 27.526  on 19  degrees of freedom
## AIC: 29.526
## 
## Number of Fisher Scoring iterations: 3
exp(coef(m))/(1 + exp(coef(m)))
## (Intercept) 
##        0.55

Intervalle de confiance

On peut également estimer un intervalle de confiance autour de cette probabilité avec la fonction confint. Différentes méthodes existent pour obtenir un intervalle de confiance à partir d’un profil de vraisemblance, mais nous n’arborderons pas cela ici (i.e. profile likelihood confidence intervals).

ci <- confint(m)  # calcul de l'intervalle de confiance sur les coefficients
## Waiting for profiling to be done...
ci <- c(coef(m), ci)  # intercept est l'estimation de la probabilité p
c(exp(ci)/(1 + exp(ci)))  # intervalle de confiance
## (Intercept)       2.5 %      97.5 % 
##   0.5500000   0.3359294   0.7519795

Log de la vraisemblance

Pour des raisons pratiques, on va souvent maximiser le log de la vraisemblance (log-likelihood) pour éviter des problèmes numériques associés aux très petites valeurs. On parle donc souvent de log-likelihood (LL, logLik, etc.).
Par exemple, on peut comparer la vraisemblance calculée par le glm à celle que nous avons calculée “à la mitaine” avec la fonction dbinom pour la valeur du paramètre qui maximise cette vraisemblance (p = 0.55).

prod(dbinom(y, size = 1, prob = 0.55))
## [1] 1.05415e-06
exp(logLik(m))
## 'log Lik.' 1.05415e-06 (df=1)

Probabilités conditionnelles

Voyons d’abord ce qu’est une probabilité conditionnelle.
\[P(A|B)\]
Ici, on représente la probabilité de \(A\) considérant que \(B\) s’est produit. Par exemple, quelle est la probabilité d’être infecté par la Covid-19 si on reçoit un test positif.

\[P(\textrm{infecté}|\textrm{test positif})\]

On pourrait aussi s’intéresser à la probabilité d’avoir un test positif si on est effectivement infecté. \[P(\textrm{test positif}|\textrm{infecté})\]

Règle de Bayes

Il s’agit simplement d’une règle en théorie des probabilités qui permet de relier des probabilités conditionnelles.

\[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\]

\(P(A)\) veut dire la probabilité de \(A\), indépendemment de \(B\).

\(P(A|B)\) veut dire la probabilité de \(A\) sachant \(B\). C’est une probabilité conditionnelle.

Interprétation

Un exemple classique de cette règle est le suivant. Supposons que \(P(A)\) désigne la probabilité d’être infecté par la Covid-19 et que \(P(B)\) désigne la probabilité de tester positif à un test de Covid-19. Ici, on a aussi \(P(A')\) qui est la probabilité de ne pas être infecté et \(P(B')\) qui est la probabilité d’un test négatif.

Si on a:
\(P(A) = 1/500 = 0.002\): proportion de personnes infectées

\(P(A') = 1-1/500 = 0.998\): proportion de personnes qui ne sont pas infectées

\(P(B|A) = 0.98\): probabilité d’un test positif si une personne est infectée

\(P(B|A') = 0.05\): probabilité d’un test positif si une personne n’est pas infectée

\(P(B) = P(B|A)P(A) + P(B|A')P(A')\): il faut additionner les possibilités associées au cas où une personne est infectée et au cas où une personne n’est pas infectée.


Quelle est la probabilité d’être infecté si on reçoit un test positif (en assumant qu’on se fait tester pour une raison autre que d’avoir des symptômes) ?

\[P(A|B) = \frac{P(B|A)P(A)}{P(B|A)P(A)+P(B|A')P(A')} = \frac{0.98 \times 0.002}{(0.98 \times 0.002) + (0.05 \times 0.998)} \approx 0.0378\]


Dans ce cas fictif, la probabilité d’avoir réellement la Covid-19 si le test est positif n’est que d’environ 3.78 %

Statistiques bayésiennes

Les statisiques bayésiennes se basent sur cette règle pour faire une inteprétation probabiliste des paramètres à partir des données.
\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]

\(data\): données
\(\theta\): paramètres du modèle

\(P(\theta|data)\): distribution postérieure (posterior)
\(P(data|\theta)\): vraisemblance (likelihood)
\(P(\theta)\): distribution a priori (prior)
\(P(data)\): distribution marginale des données (marginal likelihood)

Distribution postérieure (posterior)

\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]

La distribution postérieure \(P(\theta|data)\) est ce que l’on cherche à estimer. Il s’agit d’une distribution de probabilité qui décrit les probabilités des différentes valeurs possibles du paramètre. Contrairement aux analyses fréquentistes, on peut directement parler de la probabilité qu’un paramètre soit entre telle et telle valeur.

Vraisemblance (likelihood)

\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]

Comme vu précédemment, la vraisemblance \(P(data|\theta)\) représente la probabilité des données conditionnelles à des valeurs de paramètres \(\theta\).

Distribution a priori (priors)

\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]

Les distributions a priori \(P(\theta)\) représentent nos croyances ou nos connaissances sur les valeurs des paramètres avant de faire les analyses. On les représente à l’aide de distributions de probabilité. Ces croyances représentent ce à quoi on s’attend en termes de valeurs en se basant sur nos connaissances avant de faire notre expérience.

Distribution marginale des données

\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)} = \frac{P(data|\theta)P(\theta)}{\int P(data|\theta)P(\theta)d\theta}\]

Cette quantité \(P(data)\) (marginal likelihood) est beaucoup plus difficile à comprendre intuitivement. Elle représente en quelque sorte la probabilité des données pour l’ensemble des valeurs de paramètres possibles (d’où l’intégrale). C’est la distribution marginale des données ou la vraisemblance marginale. Elle permet également de faire en sorte que la distribution postérieure soit une distribution de probabilité (i.e. aire sous la courbe = 1).

Il s’agit d’une constante qui est généralement très difficile à calculer. Des méthodes comme le MCMC permettent d’éviter de la calculer directement dans l’obtention de la distribution postérieure.

En résumé

\[Posterior = \frac{Likelihood \times Prior}{Marginal\;Likelihood}\]



Puisque le dénominateur est une constante et que dans bien des cas, on peut éviter de le calculer, on peut plus simplement écrire:



\[Posterior \propto Likelihood \times Prior\]


En d’autres mots, la distribution postérieure est proportionnelle à la vraisemblance et aux distributions a priori



On pourrait aussi dire qu’avec une approche bayésienne, on met à jour nos croyances ou nos connaissances à l’aide de données.

Avantages

  • Interprétation probabiliste plus facile (e.g. intervalle de confiance fréquentiste vs. bayésien)

  • Facilite l’estimation de modèles complexes (e.g. modèle hiérarchique)

  • Intégration plus facile de l’incertitude à toutes les étapes du processus

  • Facilite l’estimation de quantités dérivées (e.g. prédictions ou combinaisons de prédictions)

  • Le fait d’avoir à spécifier des priors nous force à comprendre ce que nos paramètres représentent

  • Façon de penser plus naturelle, philosophiquement plus satisfaisante?

Désavantages ?

  • Souvent plus complexe (mais pas toujours)

  • Modèles plus longs à estimer

  • Moins de solutions clés en main?

  • Nécessité de choisir et de justifier ses priors

Exemple

Simulons d’abord un modèle avec des paramètres connus.

set.seed(1234)
n <- 100
x <- runif(n, 0, 10)
y <- 10 + 2 * x + rnorm(n, 0, 10)
d <- data.frame(y, x)

Quels priors utiliser?

Il est toujours préférable de visualiser les priors et de réfléchir à ce qu’ils représentent.

par(mfrow = c(1, 3))
curve(dnorm(x, 0, 50), from = -200, to = 200, main = "Intercept normal(0,50)", ylab = "Density")
curve(dnorm(x, 0, 5), from = -20, to = 20, main = "Pente normal(0,5)", ylab = "Density")
curve(dgamma(x, 2, 0.1), from = 0, to = 100, main = "Écart-type gamma(2,0.1)", ylab = "Density")


Par exemple, pour la variance/écart-type du modèle, il faut idéalement un prior avec des valeurs strictement positives, car les valeurs négatives ne font pas de sens pour la variance. Une des possibilités est d’utiliser la distribution gamma.

Modèle avec brms

La façon la plus “clé en main” de faire un modèle bayésien est d’utiliser le package brms. Ce package permet de spécifier des modèles d’une façon très similaire aux packages fréquemment utilisés comme nlme, lme4, glmmTMB, etc. En arrière-plan, brms utilise le programme STAN qui est un environnement de programmation pour la modélisation bayésienne.

library(brms)

m <- brm(y ~ x,
       data = d,
       cores = 3,
       chains = 3,
       family = "gaussian",
       prior = c(
         set_prior("normal(0,5)", class = "b", coef="x"),
         set_prior("normal(0,50)", class = "Intercept"),
         set_prior("gamma(2,0.1)", class = "sigma")
       )
     )

Visualisation

library(ggeffects)
g <- ggpredict(m, terms = c("x"))
plot(g, add.data = TRUE)

Distribution postérieure des prédictions

newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))

pp <- posterior_epred(m, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Formulation du modèle

On a la composante aléatoire, la composante systématique et les priors.


\[y \sim \mathcal{N}(\mu,\sigma^2)\] \[\mu = \beta_{0}+\beta_{1}x\] \[\beta_{0} \sim \mathcal{N}(0,50)\] \[\beta_{1} \sim \mathcal{N}(0,5)\] \[\sigma \sim \mathcal{Gamma}(2,0.1)\]

Sommaire

summary(m)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ x 
##    Data: d (Number of observations: 100) 
##   Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 3000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    10.97      1.85     7.25    14.52 1.00     2892     2159
## x             1.96      0.36     1.26     2.69 1.00     2997     2243
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     9.68      0.72     8.42    11.19 1.00     2611     1835
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).


Notez encore une fois ici qu’on peut faire une interprétation probabiliste plus naturelle des intervalles de confiance. Contrairement aux analyses fréquentistes, on peut dire que le probabilité que la valeur d’un paramètre soit entre les bornes inférieures et supérieures d’un intervalle de confiance à 95% est de 95%. En fait, avec une approche bayésienne, on parle plutôt d’intervalle de crédibilité (credible interval).

Visualisation des paramètres

plot(m)

Les distributions postérieures des paramètres sont illustrées à gauche. À droite, les chaînes de MCMC sont illustrées (nous y reviendrons).

MCMC

Le MCMC (Markov Chain Monte Carlo) est une classe d’algorithmes permettant d’échantillonner à partir d’une distribution de probabilité. Il s’agit d’une approche qui permet d’éviter des calculs complexes, voire impossible, en utilisant une approche d’échantillonnage.

Dans un contexte bayésien, le MCMC permet d’échantillonner à partir d’une distribution postérieure. Autrement dit, la distribution postérieure de chaque paramètre estimé dans un contexte bayésien peut entre autre être obtenue par cette méthode d’échantillonnage.

Sous certaines conditions, les valeurs échantillonnées représenteront la distribution postérieure de chaque paramètre.

Échantillonnage des chaînes

L’échantillonnage par MCMC se fait de façon séquentielle, i.e. que les valeurs seront échantillonnées les unes à la suite des autres. Dans bien des cas, les valeurs actuelles seront dépendentes (corrélées) des valeurs précédentes. On peut parler ici de chaînes pour désigner l’ensemble des valeurs échantillonnées de façon séquentielle.

Pour démarrer chaque chaîne, il faut d’abord sélectionner une valeur de départ pour chaque paramètre estimé. Cette valeur de départ peut potentiellement être très loin d’une valeur probable.

Éventuellement, les chaînes vont s’éloigner de ces valeurs de départ pour converger vers les distributions postérieures. À ce moment, on dira que l’échantillonnage se fait véritablement à partir des distributions postérieures.

plot(m)

Convergence des chaînes

Avec les analyses bayésiennes où les paramètres sont estimés par MCMC, il faut donc vérifier la convergence des chaînes pour s’assurer que les valeurs échantillonnées représentent bien les distributions postérieures.

On élimine souvent une certaine portion du début des chaînes pour enlever les valeurs échantillonnées, alors que les chaînes n’ont pas convergé. On parle souvent de burnin pour désigner cette partie des chaînes.

Il est également très fréquent d’utiliser plusieurs chaînes (3-4) pour un même paramètre ce qui permet de s’assurer que les chaînes ont toutes convergé au même endroit. En général, on choisira des valeurs de départ différentes pour chacune des chaînes.

On peut vérifier la convergence des chaînes avec entre autres le Rhat (potential scale reduction statistic ou Gelman-Rubin statistic). Intuitivement, pour un paramètre donné, cette valeur compare la variance inter-chaîne à la variance intra-chaîne pour évaluer si les chaînes ont convergé au même endroit. Si c’est le cas, le Rhat sera très près de 1 (Rhat <= 1.01)

summary(m)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ x 
##    Data: d (Number of observations: 100) 
##   Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 3000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    10.97      1.85     7.25    14.52 1.00     2892     2159
## x             1.96      0.36     1.26     2.69 1.00     2997     2243
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     9.68      0.72     8.42    11.19 1.00     2611     1835
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Priors

L’utilisation de priors est souvent vue comme étant la partie la plus controversée des statistiques bayésiennes. Qu’en pensez-vous?

Idéalement, il faut toujours avoir une bonne idée de l’influence des priors sur nos analyses. Pour cela il est bien de comparer les distributions postérieures aux priors. On peut également faire dans certains cas des analyses de sensibilité pour quantifier cette influence.

Plus on a de données, moins les priors auront d’influence sur l’estimation des paramètres.

Types de priors

  • Informatif: basé sur la théorie, reflète les connaissances.

  • Non-informatif: reflète l’absence de connaissances ou le désir de ne pas influencer les résultats.

  • Plat (flat): similaire à non-informatif.

  • Légèrement informatif (weaky informative): la théorie n’est pas assez développée pour utiliser des priors informatifs, mais assez pour savoir que certaines valeurs ne sont pas probables.

  • Par défaut: valeurs par défaut utilisées par les logiciels ou les packages.

  • Régularisant (regularizing): permet de circonscrire les valeurs obtenues et de pénaliser les valeurs élevées dans un contexte prédictif.

  • Etc.


Voir notamment Lemoine 2019 et Banner et al. 2020 pour une discussion sur l’utilisation des priors en écologie.

Quels priors utiliser?

Il faut aussi sélectionner une distribution approriée pour représenter nos priors et les distributions utilisées doivent être appropriées pour un paramètre donné. Le groupe qui travaille sur la plateforme STAN a plusieurs recommandations de priors à utiliser dépendemment du contexte (priors pour des pentes, priors pour des variances, priors pour des variances dans un contexte de modèles hiérarchiques, etc.)

Voir Prior Choice Recommendations pour des recommandations.

Attention aux priors “non-informatifs”!

Les priors non-informatifs ne le sont pas toujours sous l’échelle du paramètre. Par exemple, si vous tentez d’estimer une probabilité et que vous utiliser une distribution normale comme prior pour l’intercept du modèle, il se peut que ce prior soit très informatif sous l’échelle du prédicteur linéaire qui est relié à la probabilité par une fonction logit. Supposons un glm binomial avec seulement un intercept.

\[ y \sim \mathcal{Binom}(1,p) = \mathcal{Bernoulli}(p)\]
\[ logit(p) = log\left(\frac{p}{1-p}\right) = \beta_{0} \]
\[ \beta_{0} \sim \mathcal{N}(0,10) \]

Attention aux priors “non-informatifs”!

n <- 10000
x <- rnorm(n, 0, 10)
par(mfrow = c(1, 2))
hist(x, breaks = 50, main = "Valeurs a priori de l'intercept")
hist(exp(x)/(1 + exp(x)), breaks = 50, main = "Lorsque converties en probabilité")


Autrement dit, le prior est très informatif et pousse les valeurs vers des probabilités extrêmes!

Étudier ses priors

La distribution prédictive préalable (Prior predictive distribution) permet de visualiser l’effet des priors sur l’échelle de la réponse.

m2 <- update(m, sample_prior = "only")
newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))

pp <- posterior_epred(m2, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Avec un prior aberrant

newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))

m3 <- update(m, 
        sample_prior = "only",
        prior = c(
           set_prior("normal(50,20)", class = "b", coef="x"),
           set_prior("normal(0,50)", class = "Intercept"),
           set_prior("gamma(2,0.1)", class = "sigma")
        ),
      )
pp <- posterior_epred(m3, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Avec un modèle binomial

logit <- function(p) {
  log(p / (1 - p))
}

inv.logit <- function(x) {
   exp(x) / (1 + exp(x))  
}

n <- 200
x <- runif(n, 0, 10)
z <- -2 + 0.5 * x
y <- rbinom(n, size = 1, prob = inv.logit(z))

d<-data.frame(y,x)

m <- brm(y ~ x,
       data = d,
       cores = 3,
       chains = 3,
       family = "bernoulli",
       prior = c(
         set_prior("normal(0,5)", class = "b", coef="x"),
         set_prior("normal(0,5)", class = "Intercept")
       )
     )

Distribution postérieure des prédictions

newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))

pp <- posterior_epred(m, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Ce qu’indique notre prior

m2 <- update(m, sample_prior = "only")
pp <- posterior_epred(m2, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Prédiction pour x = 2

pp <- posterior_epred(m2, data.frame(x = 2), ndraws = 500)
hist(pp, breaks = 50)

Un prior plus raisonnable?

m3 <- update(m,
       sample_prior = "only",       
       prior = c(
         set_prior("normal(0,1)", class = "b", coef="x"),
         set_prior("normal(0,1)", class = "Intercept")
       )
     )
pp <- posterior_epred(m3, newdata, ndraws = 500)

plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35)) 
lapply(1:nrow(pp), function(i){
  lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()

Outils

Il existe plusieurs outils pour faire des modèles bayésiens.

WinBUGS Un des premiers outils, mais un peu lent

JAGS Similaire à BUGS, mais un peu plus rapide

STAN Rapide et bien documenté

Nimble Relativement nouveau

MCMCglmm Modèle animal?

brms Facile d’utilisation et basé sur STAN

INLA Modèles spatiaux